///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_REDUCTIONS_H_
#define _BUOLA_MAT_REDUCTIONS_H_

#include <buola/mat/cexpr.h>
#include <buola/mat/detail/reductions.h>

namespace buola { namespace mat {

namespace detail {

template<typename tReduction,typename tExpr>
OScalar<tExpr> deduct_reduction_result(const CExpr<tExpr>&);

template<typename tReduction,typename tSub>
detail::CColReduction<const tSub&,tReduction> deduct_reduction_result(const CColwiseOp<tSub>&);

template<typename tReduction,typename tSub>
detail::CRowReduction<const tSub&,tReduction> deduct_reduction_result(const CRowwiseOp<tSub>&);

template<typename tExpr,typename tReduction>
using OReductionResult=decltype(deduct_reduction_result<tReduction>(std::declval<const tExpr&>()));

template<typename tReduction,typename tScalar,typename tExpr>
tScalar deduct_reduction_result_with_type(const CExpr<tExpr>&);

template<typename tReduction,typename tScalar,typename tSub>
detail::CColReduction<const tSub&,tReduction,tScalar> deduct_reduction_result_with_type(const CColwiseOp<tSub>&);

template<typename tReduction,typename tScalar,typename tSub>
detail::CRowReduction<const tSub&,tReduction,tScalar> deduct_reduction_result_with_type(const CRowwiseOp<tSub>&);

template<typename tExpr,typename tReduction,typename tScalar>
using OReductionResultWithType=decltype(deduct_reduction_result_with_type<tReduction,tScalar>(std::declval<const tExpr&>()));

template<typename tExpr,typename tReduction>
detail::OScalar<tExpr> do_reduction(const CExpr<tExpr> &pE,const tReduction &pR)
{
    return pR(detail::element_accessible<detail::ACCESS_ONCE>(pE.Derived()));
}

template<typename tSub,typename tReduction>
std::pair<const tSub&,tReduction> do_reduction(const detail::CColwiseOp<tSub> &pE,const tReduction &pR)
{
    return {pE.Sub(),pR};
}

template<typename tSub,typename tReduction>
std::pair<const tSub&,tReduction> do_reduction(const detail::CRowwiseOp<tSub> &pE,const tReduction &pR)
{
    return {pE.Sub(),pR};
}

template<typename tScalar,typename tExpr,typename tReduction>
tScalar do_reduction_with_type(const CExpr<tExpr> &pE,const tReduction &pR)
{
    return pR(detail::element_accessible<detail::ACCESS_ONCE>(pE.Derived()));
}

template<typename tScalar,typename tSub,typename tReduction>
std::pair<const tSub&,tReduction> do_reduction_with_type(const detail::CColwiseOp<tSub> &pE,const tReduction &pR)
{
    return {pE.Sub(),pR};
}

template<typename tScalar,typename tSub,typename tReduction>
std::pair<const tSub&,tReduction> do_reduction_with_type(const detail::CRowwiseOp<tSub> &pE,const tReduction &pR)
{
    return {pE.Sub(),pR};
}

/*namespace detail*/ }
    
template<typename tDerived>
detail::CColwiseOp<const tDerived&> cols(const CExpr<tDerived> &pExpr)
{
    return {pExpr.Derived()};
}
    
template<typename tDerived>
detail::CRowwiseOp<const tDerived&> rows(const CExpr<tDerived> &pExpr)
{
    return {pExpr.Derived()};
}
    
template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionSum> sum(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionSum())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionProd> prod(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionProd())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionMean> mean(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionMean())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionNorm> norm(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionNorm())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionNorm2> norm2(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionNorm2())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionMin> min(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionMin())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionMax> max(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionMax())};
}

template<typename tExpr>
detail::OReductionResultWithType<tExpr,detail::CReductionAny,bool> any(const tExpr &pE)
{
    return {detail::do_reduction_with_type<bool>(pE,detail::CReductionAny())};
}

template<typename tExpr>
detail::OReductionResultWithType<tExpr,detail::CReductionAll,bool> all(const tExpr &pE)
{
    return {detail::do_reduction_with_type<bool>(pE,detail::CReductionAll())};
}

template<typename tExpr>
detail::OReductionResultWithType<tExpr,detail::CReductionCount,size_t> count(const tExpr &pE)
{
    return {detail::do_reduction_with_type<size_t>(pE,detail::CReductionCount())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionVar> var(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionVar())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionStdDev> stddev(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionStdDev())};
}

template<typename tExpr>
detail::OReductionResult<tExpr,detail::CReductionLogSumExp> logsumexp(const tExpr &pE)
{
    return {detail::do_reduction(pE,detail::CReductionLogSumExp())};
}

template<typename tDerived,typename tPredicate>
detail::OScalar<tDerived> sum_if(const CExpr<tDerived> &pE,tPredicate &&pPred)
{
    ///\todo optimize this
    const auto &lE=detail::element_accessible<detail::ACCESS_ONCE>(pE.Derived());
    
    detail::OScalar<tDerived> lSum=0.0;
    
    for(int r=0;r<lE.Rows();r++)
    {
        for(int c=0;c<lE.Cols();c++)
        {
            if(pPred(lE(r,c)))
                lSum+=lE(r,c);
        }
    }    
    
    return lSum;
}

template<typename tDerived>
detail::OScalar<tDerived> trace(const CExpr<tDerived> &pE)
{
    const auto &lE=detail::element_accessible<detail::ACCESS_ONCE>(pE.Derived());
    int lDim=std::min(lE.Rows(),lE.Cols());

    detail::OScalar<tDerived> lSum=0.0;
    
    for(int i=0;i<lDim;i++)
        lSum+=lE(i,i);

    return lSum;
}

template<typename tExpr>
detail::OReductionResultWithType<tExpr,detail::CReductionIndMax,size_t> indmax(const tExpr &pE)
{
    return {detail::do_reduction_with_type<size_t>(pE,detail::CReductionIndMax())};
}

template<typename tExpr>
detail::OReductionResultWithType<tExpr,detail::CReductionIndMin,size_t> indmin(const tExpr &pE)
{
    return {detail::do_reduction_with_type<size_t>(pE,detail::CReductionIndMin())};
}

template<typename tLH,typename tRH>
detail::OScalar<tLH> distance(const CExpr<tLH> &pLH,const CExpr<tRH> &pRH)
{
    return norm(pRH-pLH);
}

template<typename tLH,typename tRH>
detail::OScalar<tLH> distance2(const CExpr<tLH> &pLH,const CExpr<tRH> &pRH)
{
    return norm2(pRH-pLH);
}

template<typename tDerived>
detail::OMatrixType<tDerived> normalize(const CExpr<tDerived> &pE)
{
    return pE/norm(pE);
}

template<typename tDerived,typename tPredicate>
size_t count_if(const CExpr<tDerived> &pE,tPredicate &&pPred)
{
    const auto &lE=detail::element_accessible<detail::ACCESS_ONCE>(pE.Derived());
    
    size_t lCount=0;
    
    for(int r=0;r<lE.Rows();r++)
    {
        for(int c=0;c<lE.Cols();c++)
        {
            if(pPred(lE(r,c)))
                lCount++;
        }
    }    
    
    return lCount;
}

/*namespace mat*/ } /*namespace buola*/ }

#endif
